R Notebook

library(ggplot2)

ll.nbinom <- function(par, x, d) {
  if (min(par) < 0) {
    return(Inf)
  }
  -sum(dnbinom(x, par[1], (par[2]/d)/((par[2]/d) + 1), log=TRUE))
}

lk.nbinom <- function(par, x, d, priors=list(a=c(1, 0.5), b=c(1, 0.5))) {
  if (min(par) < 0) {
    return(-Inf)
  }
  r <- par[1]
  p <- par[2]/(par[2] + d)
  sum(dnbinom(x, r, p, log=TRUE)) + dgamma(par[1], priors$a[1], priors$a[2], log=TRUE) + dgamma(par[2], priors$b[1], priors$b[2], log=TRUE)
}

rtnorm <- function(n, mean = 0, sd = 1, lower=-Inf, upper=Inf) {
  # range is a vector of two values
  
  F.a <- pnorm(lower, mean = mean, sd = sd)
  F.b <- pnorm(upper, mean = mean, sd = sd)
  
  u <- runif(n, min = F.a, max = F.b)
  
  qnorm(u, mean = mean, sd = sd)
}

metropolis <- function(target, initial, iterations=10000, warmup=2000, thin=100) {
  k <- length(initial)
  chain <- matrix(0, iterations + warmup, k)
  accepted <- 0
  
  chain[1, ] <- initial
  
  for (i in 2:(iterations + warmup)) {
    current <- chain[i - 1, ]
    candidate <- rnorm(k, mean=current, sd=1)
    # candidate <- rtnorm(k, mean=current, sd=1, lower=0, upper=)
    A <- exp(target(candidate) - target(current))
    if (runif(1) < A) {
      accepted <- accepted + 1
      chain[i, ] <- candidate
    } else {
      chain[i, ] <- current
    }
  }
  # print(accepted/(iterations+warmup))
  return(chain[seq(warmup, iterations + warmup, thin), ])
}

gibbs <- function(target, initial, iterations=10000, warmup=2000, thin=100) {
  k <- length(initial)
  chain <- matrix(0, iterations + warmup, k)
  accepted <- 0
  
  chain[1, ] <- initial
  
  for (i in 2:(iterations + warmup)) {
    current <- chain[i - 1, ]
    candidate <- chain[i - 1, ]
    for (j in 1:k) {
      candidate[j] <- rnorm(1, mean=current[j], sd=1)
      # candidate[j] <- rtnorm(1, mean=current[j], sd=1, lower=0, upper=Inf)
      A <- exp(target(candidate) - target(current))
      if (runif(1) < A) {
        current <- candidate
      } else {
        candidate <- current
      }
    }
    chain[i, ] <- current
  }
  # print(accepted/(iterations+warmup))
  return(chain[seq(warmup, iterations + warmup, thin), ])
}

plot(density(rgamma(1000, 1, 1)))
lines(density(rgamma(1000, 1, 0.1)), col=2)

n <- 100
iterations <- 100
cors_true <- rep(0, iterations)
cors_est <- rep(0, iterations)
cors_est_lb <- rep(0, iterations)
cors_est_ub <- rep(0, iterations)
mean_sampling <- rep(0, iterations)
bad_chains <- rep(0, iterations)

i <- 1

for (i in 1:iterations) {
  print(i)
  a <- runif(1, 1, 10)
  b <- runif(1, 1, 10)
  
  # Do variable sampling effort.
  mean_sampling[i] <- runif(1, 1, 20)
  d <- rpois(n, mean_sampling[i])
  d[d == 0] = 1
  
  alpha <- rgamma(n, a, b)
  alpha_hat <- rpois(n, alpha*d)/d
  
  cors_true[i] <- cor(alpha, alpha_hat)
  # cors_true[i] <- var(alpha)/var(alpha_hat)

  # parameters <- optim(c(1, 1), ll.nbinom, x=as.integer(alpha_hat * d), d=d)$par
  # a_ <- parameters[1]
  # b_ <- parameters[2]
  
  target <- function(par) lk.nbinom(par, as.integer(alpha_hat * d), d)
  chain <- gibbs(target, c(1, 1), warmup=2000, iterations=20000, thin=10)

  plot(chain[, 1], type="l", col=1, ylim=c(0, max(chain)))
  lines(chain[, 2], type="l", col=2)
    
  qt <- quantile(chain[, 2], probs=c(0.025, 0.5, 0.975))

  b_ <- qt[2]
  b_lb <- qt[1]
  b_ub <- qt[3]
  
  cors_est[i] <- sqrt(n/(n + b_ * sum(d^-1)))
  cors_est_ub[i] <- sqrt(n/(n + b_ub * sum(d^-1)))
  cors_est_lb[i] <- sqrt(n/(n + b_lb * sum(d^-1)))
}

[1] 1
[1] 2

[1] 3

[1] 4

[1] 5

[1] 6

[1] 7

[1] 8

[1] 9

[1] 10

[1] 11

[1] 12

[1] 13

[1] 14

[1] 15

[1] 16

[1] 17

[1] 18

[1] 19

[1] 20

[1] 21

[1] 22

[1] 23

[1] 24

[1] 25

[1] 26

[1] 27

[1] 28

[1] 29

[1] 30

[1] 31

[1] 32

[1] 33

[1] 34

[1] 35

[1] 36

[1] 37

[1] 38

[1] 39

[1] 40

[1] 41

[1] 42

[1] 43

[1] 44

[1] 45

[1] 46

[1] 47

[1] 48

[1] 49

[1] 50

[1] 51

[1] 52

[1] 53

[1] 54

[1] 55

[1] 56

[1] 57

[1] 58

[1] 59

[1] 60

[1] 61

[1] 62

[1] 63

[1] 64

[1] 65

[1] 66

[1] 67

[1] 68

[1] 69

[1] 70

[1] 71

[1] 72

[1] 73

[1] 74

[1] 75

[1] 76

[1] 77

[1] 78

[1] 79

[1] 80

[1] 81

[1] 82

[1] 83

[1] 84

[1] 85

[1] 86

[1] 87

[1] 88

[1] 89

[1] 90

[1] 91

[1] 92

[1] 93

[1] 94

[1] 95

[1] 96

[1] 97

[1] 98

[1] 99

[1] 100

qplot(cors_true, cors_est) +
  geom_errorbar(aes(ymin=cors_est_lb, ymax=cors_est_ub)) +
  geom_smooth() +
  geom_abline() +
  xlab("True correlation") +
  ylab("Estimated correlation")

qplot(mean_sampling, abs(cors_true - cors_est)) +
  # geom_errorbar(aes(ymin=abs(cors_true-cors_est_lb), ymax=abs(cors_true-cors_est_ub))) +
  geom_smooth()

a <- runif(1, 1, 10)
b <- runif(1, 1, 10)

X <- seq(0.1, 100, 0.1)

alpha <- rgamma(n, a, b)

priors <- list(a=c(1, 1), b=c(1, 1))

d <- rep(1, n)
alpha_hat <- rpois(n, alpha*d)/d
plot(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=1)

d <- rep(10, n)
alpha_hat <- rpois(n, alpha*d)/d
lines(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=2)

d <- rep(100, n)
alpha_hat <- rpois(n, alpha*d)/d
lines(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=3)